#include <cmath>
#include <iostream>
#include <time.h>
#include "randomc.h"
#include "MT.hpp"
#include "atom.h"
#include <fstream>
#define PI 3.1415926535897932384626433832795

int nint(double x)
{
	int i = int(x);
	if (double(x-i)>0.5)
		return i+1;
	else if (double(x-i)<-0.5)
		return i-1;
	else return i;
}

Atom::Atom(const unsigned& n, const unsigned &m)
{
	N = n;
	Nbond=m;
	mass = new double[N];
	x = new double*[N];
	v = new double*[N];
	xrel = new double*[N];
	vrel = new double*[N];
	f = new double*[N];
	boxfold = new int*[N];

	for(unsigned i = 0; i < N; ++i)
	{
		mass[i] = 1.0;
		x[i] = new double[3];
		v[i] = new double[3];
		f[i] = new double[3];
		xrel[i] = new double[3];
		vrel[i] = new double[3];
		boxfold[i] = new int[3];
	}

	for(unsigned i = 0; i < N; ++i)
	for(unsigned j = 0; j < 3; ++j)
		boxfold[i][j] = 0;

	Bond_ends = new unsigned*[Nbond];
	Bond_length = new double[Nbond];
			
	for(unsigned i = 0; i < Nbond; ++i)
		Bond_ends[i]=new unsigned[2];
}

Atom::~Atom()
{
	delete mass;
	delete x;
	delete v;
	delete f;
	delete group;
	delete Mass;
	delete boxfold;
	delete Bond_ends;
	delete Bond_length;
	delete xcen;
	delete vcen;
	delete fcen;
	delete xrel;
	delete vrel;
	delete diff;
	delete xcen0;
	delete r_dis;
	delete x0;
	delete []vcen0;
	delete []v_corr;
}

void Atom::group_init()
{

	Nm = new unsigned[Ngroup];
	Nstart=new unsigned[Ngroup];
	rg = new double[Ngroup];
	group = new unsigned*[Ngroup];
	Mass = new double[Ngroup];
	xcen = new double*[Ngroup];
	vcen = new double*[Ngroup];
	fcen = new double*[Ngroup];

	for(unsigned i = 0; i < Ngroup; ++i)
	{
		if(i<Ngroup_p)
		{
			Nm[i]=Nm_p;
			rg[i] = rg_p;
			Nstart[i]=i*Nm_p;
		}
		else
		{
			Nm[i]=Nm_s;
			rg[i] = rg_s;
			Nstart[i]=Ngroup_p*Nm_p+(i-Ngroup_p)*Nm_s;
		}
	  
		group[i] = new unsigned[Nm[i]];
		xcen[i] = new double[3];
		vcen[i] = new double[3];
		fcen[i] = new double[3];
		Mass[i] = 0.0;
    }
 
	for(unsigned i = 0; i < Ngroup; ++i)
	for(unsigned j = 0; j < Nm[i]; ++j)
	{
		group[i][j] = Nstart[i]+j;
		Mass[i] += *(mass + Nstart[i]+j );
	}
	
	if(Nstart[Ngroup-1]+Nm[Ngroup-1] != N)
	{
		std::cout<<" k=\t"<<Nstart[Ngroup-1]+Nm[Ngroup-1]<<"\n N=\t"<<N;
		std::cout<<"\ncount number error, please check! @atom.cpp L115.\n";
		std::cin.get();
	}

	std::cout<<"Matrix defination finished. \n";
  
  //diffusion viarable
	unsigned nd = dft/dfdt;
	std::cout<<"nd=dft/dfdt=\t"<<nd<<"\t@atom.cpp L130.\n";
	diff = new double*[nd];
	visc = new double*[nd];
	v_corr = new double*[nd];
  
	for(unsigned i = 0; i <nd; ++i)
	{
		diff[i] = new double[2];
		visc[i] = new double[2];
	}
  
	for(unsigned i = 0; i < 10 * nd; ++i)
		v_corr[i] = new double[2];				//  ???
 
	for(unsigned i = 0; i < nd; ++i)
	for(unsigned j = 0; j < 2; ++j)
	{
		diff[i][j] = 0;
		visc[i][j] = 0;
	}
  
	xcen0 = new double*[Ngroup];
	vcen0 = new double*[Ngroup];
	for(unsigned i = 0 ;i < Ngroup; ++i)
	{
		xcen0[i] = new double[3];
		vcen0[i] = new double[3];
	}
  
	//neighbor particle distribution 
	particle_number = new unsigned[Ngroup];
	pnumber_dis = new unsigned[Ngroup];
  
	x0 = new double*[N];
	for(unsigned i = 0; i < N; ++i)
		x0[i] = new double[3];
  
	r_dis = new double[N_radial];
	for(int i = 0; i < N_radial; ++i)
		r_dis[i] = 0;
	r_dis_num = 0;

	veldis = new double*[velocitybin];
	for(int i = 0; i < velocitybin; ++i)
		veldis[i] = new double[2];

	for(int i = 0; i < velocitybin; ++i)
	for(int j = 0; j < 2; ++j)
		veldis[i][j] = 0;
  
	for(unsigned i = 0; i < Ngroup; ++i)
	{
		particle_number[i] = 0;
		pnumber_dis[i] = 0;
	}
}


double Atom::dis_sq(int i, int j)
{
	double xij = x[i][0] - x[j][0];
	double yij = x[i][1] - x[j][1];
	double zij = x[i][2] - x[j][2];

	xij = xij - nint(xij/boxx) * boxx;
	yij = yij - nint(yij/boxy) * boxy;
	zij = zij - nint(zij/boxz) * boxz;

	double rsq = xij*xij + yij*yij + zij*zij;
	return rsq;
}

double Atom::temper()
{
	double sum = 0;
	double sumv[3] = {0,0,0};
	
	for(unsigned i = 0; i < N; ++i)
	for(unsigned j = 0; j < 3; ++j)
		sumv[j] += v[i][j] / N;

	for(unsigned i = 0; i < N; ++i)
	for(unsigned j = 0; j < 3; ++j)  //{
		sum += (v[i][j] - sumv[j]) * (v[i][j] - sumv[j]);
	
	return sum/(3.0*N - Ngroup);
}

void Atom::period_bd()
{
	for(unsigned i = 0; i < N; ++i)
	{
		boxfold[i][0] = nint(x[i][0]/boxx);
		boxfold[i][1] = nint(x[i][1]/boxy);
		boxfold[i][2] = nint(x[i][2]/boxz);
		x[i][0] = x[i][0] - nint(x[i][0]/boxx) * boxx;
		x[i][1] = x[i][1] - nint(x[i][1]/boxy) * boxy;
		x[i][2] = x[i][2] - nint(x[i][2]/boxz) * boxz;
	}
}


void Atom::x_init()
{
	int32 seed = (int32) time(0);
	TRandomMersenne RG(seed);
	for(unsigned i = 0; i < N; ++i)
	{
		x[i][0] = (RG.Random() - 0.5) * boxx;
		x[i][1] = (RG.Random() - 0.5) * boxy;
		x[i][2] = (RG.Random() - 0.5) * boxz;
	}  
}


void Atom::x_group_init()
{
	int32 seed =(int32) time(0);   
	TRandomMersenne RG(seed);
	MT::init_genrand(seed);

	double **x0;
	x0 = new double*[Ngroup];
	for(unsigned i = 0; i < Ngroup; ++i)
		x0[i] = new double[3];

	double r, angle1, angle2;

	if(Ngroup_p>0)
	{		
		double **pos;
		pos=new double *[Nm_p];
		for(unsigned i=0; i<Nm_p; ++i)
			pos[i]=new double[3];

		unsigned **temp_ends;
		temp_ends=new unsigned *[Nm_p-1];
		for(unsigned i=0; i<Nm_p-1; ++i)
			temp_ends[i]=new unsigned[2];

		double *temp_length;
		temp_length=new double[Nm_p-1];

		unsigned k=0;
		unsigned k_bond=0;

		unsigned N_arm=12;
		unsigned end1, endC;
		unsigned kk=0;

		k=0;
		pos[0][0]=0.0;
		pos[0][1]=0.0;
		pos[0][2]=0.0;

		endC=k;

		for(unsigned i = 0; i < N_arm; ++i)
		{
			end1=endC;	

			double r0=1.615;
			double dsita=0.85*PI/N_arm;
			double sita0=i*2*PI/N_arm;

			for(unsigned j=0; j< (Nm_p-1)/N_arm; ++j)
			{
				++k;
				double sita=sita0+j*dsita;
				if(j != 0)	r0=1.147;

				do
				{
					angle1 =   MT::genrand_real2() * PI;

					pos[k][0] = pos[end1][0]  + r0 *sin(angle1) * cos(sita);
					pos[k][1] = pos[end1][1]  + r0 *sin(angle1) * sin(sita);
					pos[k][2] = pos[end1][2]  + r0 *cos(angle1);
				}while(check(k,pos, 0.66* 0.66) == 0);

				temp_ends[kk][0]=end1;
				temp_ends[kk][1]=k;
				temp_length[kk]=r0;

				++kk;
				end1=k;
			}
		}
/*	 	
	std::ofstream pos1_out("positon1.dat");
	for(unsigned k = 0; k < Nm_p; ++k)
			pos1_out<<pos[k][0]<<"\t"<<pos[k][1]<<"\t"<<pos[k][2]<<"\n";

			std::cin.get();
*/

//		for(unsigned i=0; i<Nm_p-1; ++i)
//			std::cout<<i<<"\t"<<temp_ends[i][0]<<"\t"<<temp_ends[i][1]<<"\t"<<temp_length[i]<<"\t@atom.cpp L351.\n";

//---------------------------------------------------------------------------------------------------------------
		k=0;
		unsigned Nxy=4;				// The number of polymer in one x-y section.
		double deltaZ=boxz/(Ngroup_p/Nxy);
		if(deltaZ < 1 ){ std::cout<<"deltaZ is too small. @ atom.cpp L388\n"; std::cin.get();}
		
		x0[k][0] = 0;
		x0[k][1] = 0;
		x0[k][2] = - 0.5 * boxz + 0.5*deltaZ;
		std::cout<<"ploymer group\t"<<k<<"\tis generated. @atom.cpp L338\n";

		++k;
		for(unsigned i=1; i<Nxy; ++i)
		{
			do
			{
				x0[k][0] = (MT::genrand_real2() - 0.5) * (boxx);
				x0[k][1] = (MT::genrand_real2() - 0.5) * (boxy);
				x0[k][2] = - 0.5 * boxz + 0.5*deltaZ;
			} while (check(k, x0, 2*rg_p*rg_p) == 0);
			std::cout<<"ploymer group\t"<<k<<"\tis generated. @atom.cpp L347\n";
			++k;
		}

		for(unsigned i=1; i<Ngroup_p/Nxy; ++i)
		for(unsigned j=0; j<Nxy; ++j)
		{

			x0[k][0] = x0[j][0];
			x0[k][1] = x0[j][1];
			x0[k][2] = - 0.5 * boxz + ( 0.5+i )*deltaZ;

			std::cout<<"ploymer group\t"<<k<<"\tis generated. @atom.cpp L359\n";
			++k;
		}
		if(k != Ngroup_p) {std::cout<<"The number of polymer group is wrong @atom.cpp L364\n"<<k<<"\t"<<Ngroup_p; std::cin.get(); }

//------------------------------------------------------------------------------------------------------------------		

		for(unsigned i = 0; i <Ngroup_p; ++i)
		{
			for(unsigned j=0; j<Nm[i]; ++j)
			{
				k=group[i][j];
				kk=k-*group[i];
//				std::cout<<i<<"\t"<<kk<<"\n";
				x[k][0]=x0[i][0]+pos[kk][0];
				x[k][1]=x0[i][1]+pos[kk][1];
				x[k][2]=x0[i][2]+pos[kk][2];
			}

			for(unsigned j=0; j<Nm[i]-1; ++j)
			{
				k=i*(Nm_p-1);
				Bond_ends[k+j][0]=*group[i]+temp_ends[j][0];	
				Bond_ends[k+j][1]=*group[i]+temp_ends[j][1];;
				Bond_length[k+j]=temp_length[j];
			}
		}

//		for(unsigned i=0; i<Nbond; ++i)
//			std::cout<<i<<"\t"<<Bond_ends[i][0]<<"\t"<<Bond_ends[i][1]<<"\t"<<Bond_length[i]<<"\t@atom.cpp L351.\n";
//		std::cin.get();
		
		for(unsigned i = Ngroup_p; i <Ngroup; ++i)
		{
			do
			{
				x0[i][0] = (MT::genrand_real2() - 0.5) * (boxx);
				x0[i][1] = (MT::genrand_real2() - 0.5) * (boxy);
				x0[i][2] = (MT::genrand_real2() - 0.5) * (boxz);
			} while (check(i, x0, 1.8 * rg_s * rg_s) == 0 && check2(x0[i], N_p, x,1.8*rg[i]*rg[i]) );
	
			std::cout<<"solvent group\t"<<i<<"\tis generated. @atom.cpp L401\n";
		}

	}
	else
	{
		x0[0][0] = (MT::genrand_real2() - 0.5) * (boxx);
		x0[0][1] = (MT::genrand_real2() - 0.5) * (boxy);
		x0[0][2] = (MT::genrand_real2() - 0.5) * (boxz);
		std::cout<<"solvent group\t"<<0<<"\tis generated. @atom.cpp L410\n";
		
		for(unsigned i = 1; i <Ngroup; ++i)
		{
			do
			{
				x0[i][0] = (MT::genrand_real2() - 0.5) * (boxx);
				x0[i][1] = (MT::genrand_real2() - 0.5) * (boxy);
				x0[i][2] = (MT::genrand_real2() - 0.5) * (boxz);
			} while (check(i, x0, 1.8 * rg[i] * rg[i]) == 0);
			std::cout<<"solvent group\t"<<i<<"\tis generated. @atom.cpp L420\n";
		}
	}
	
//	std::cin.get();

	unsigned k;
	for(unsigned i = Ngroup_p; i < Ngroup; ++i)
	{
		for(unsigned j = 0; j < Nm[i]; ++j)
		{
			k=group[i][j];
			do
			{
				r = rg[i] + (MT::genrand_real2() - 0.5) * 0.2;
				angle1 =   MT::genrand_real2()  * PI;
				angle2 =   MT::genrand_real2() * 2.0 * PI;

				x[k][0] = x0[i][0] + r * sin(angle1) * cos(angle2);
				x[k][1] = x0[i][1] + r * sin(angle1) * sin(angle2);
				x[k][2] = x0[i][2] + r * cos(angle1);
			} while(check(k,x, 0.66 * 0.66) == 0);
		}
		std::cout<<i<<"\t"<<k<<"\t @atom.cpp L293.\n";
	}
	
	double rg_temp;
	double r0[3];

	for(unsigned i = 0; i < Ngroup; ++i)
	{
		rg_temp = gyration(i);
		if(rg_temp==0) {std::cout<<"Error: gyration is zero @atom.cpp L454\n"<<i; std::cin.get(); }
		
		center(r0, i);

		for(unsigned j = 0; j < Nm[i]; ++j)
		{
			k=Nstart[i]+j;
			for(unsigned l = 0; l < 3; ++l)
			{
				x[k][l] = rg[i] / rg_temp * (x[k][l] - r0[l]) + r0[l];
				xrel[k][l] = x[k][l] - r0[l];
			}
		}
		
		for(unsigned l = 0; l < 3; ++l)
			xcen[i][l] = r0[l];
	}

	period_bd();
	delete x0;
}


bool Atom::check(unsigned i, double **x0, double rc)
{
	double r2;
	double rij[3];
	bool flag = true;
	
	for(unsigned l = 0; l < i; ++l)
	{
		rij[0] = x0[l][0] - x0[i][0];
		rij[0] = rij[0] - nint(rij[0] / boxx) * boxx;
		rij[1] = x0[l][1] - x0[i][1];
		rij[1] = rij[1] - nint(rij[1] / boxy) * boxy;
		rij[2] = x0[l][2] - x0[i][2];
		rij[2] = rij[2] - nint(rij[2] / boxz) * boxz;
		r2 = rij[0] * rij[0] + rij[1] * rij[1] + rij[2] * rij[2];

		if(r2 < rc)
		{
			flag = false;
			break;
		}
	}
	return flag;
}


bool Atom::check2(double *ch, unsigned i, double **x0, double rc)
{
	double r2;
	double rij[3];
	bool flag = true;
	
	for(unsigned l = 0; l < i; ++l)
	{
		rij[0] = x0[l][0] - ch[0];
		rij[0] = rij[0] - nint(rij[0] / boxx) * boxx;
		rij[1] = x0[l][1] - ch[1];
		rij[1] = rij[1] - nint(rij[1] / boxy) * boxy;
		rij[2] = x0[l][2] - ch[2];
		rij[2] = rij[2] - nint(rij[2] / boxz) * boxz;
		r2 = rij[0] * rij[0] + rij[1] * rij[1] + rij[2] * rij[2];

		if(r2 < rc)
		{
			flag = false;
			break;
		}
	}
	return flag;
}


/*
void Atom::v_init(double &Temp)
{
  double temp[3] = {0.0, 0.0, 0.0};
  for(unsigned i = 0; i < N; ++i)
    for(unsigned j = 0; j < 3; ++j)
      {
        v[i][j] = MT::genrand_real2() - 0.5;
        temp[j] += v[i][j] / N;
      }
  double sumv2[3] = {0.0,0.0,0.0};
  
  for(unsigned i = 0; i < N; ++i)
    for(unsigned j = 0; j < 3; ++j)
      {
        v[i][j] -= temp[j];
        sumv2[j] += v[i][j] * v[i][j];
      }
  for(unsigned i = 0; i < 3; ++i)
    sumv2[i] = sqrt(N * Temp / sumv2[i]);
  for(unsigned i = 0; i < N; ++i)
    for(unsigned j = 0; j < 3; ++j)
      v[i][j] *= sumv2[j];
    
}
*/


void Atom::v_group_init(double &Temp)
{
	unsigned k;
	unsigned Nm_temp=std::max(Nm_s,Nm_p);
	double **tempx=new double*[Nm_temp];
	
	for(unsigned i = 0; i < Nm_temp; ++i)
		tempx[i]=new double[3];	

	double temp[3]; 
	int32 seed = (int32) time(0); 
	//int32 seed = 10;
	TRandomMersenne RG(seed);
  
  //construct the velocity of center of mass
  
  double sumv2[3] = {0.0, 0.0, 0.0};
  double sumv[3] = {0,0,0};

	for(unsigned i = 0; i < Ngroup; ++i)
	for(unsigned j = 0; j < 3; ++j)
	{
		vcen[i][j] = 2.0*(RG.Random() - 0.5);
		sumv[j] += vcen[i][j] / Ngroup;
	}
  
	for(unsigned i = 0; i < Ngroup; ++i)
	for(unsigned j = 0; j < 3; ++j)
	{
		vcen[i][j] -= sumv[j];
		sumv2[j] += Mass[i] * vcen[i][j] * vcen[i][j];
	}
  
	for(unsigned i = 0; i < 3; ++i)
		sumv2[i] = sqrt(Ngroup * Temp / sumv2[i]);
 
	for(unsigned i = 0; i < Ngroup; ++i)
	for(unsigned j = 0; j < 3; ++j)
		vcen[i][j] *= sumv2[j];
  
  //construct the relative velocity
	for(unsigned i = 0; i < 3; ++i)
		sumv2[i] = 0;
  
	for(unsigned i = 0; i < Ngroup; ++i)
	{
		for(unsigned j = 0; j < Nm[i]; ++j)
		{
			k=group[i][j];
			tempx[j][0] = x[k][0] + boxfold[k][0] * boxx;
			tempx[j][1] = x[k][1] + boxfold[k][1] * boxy;
			tempx[j][2] = x[k][2] + boxfold[k][2] * boxz;
		}

		for(unsigned j = 0; j < 3; ++j)
			temp[j] = 0.0;

		for(unsigned j = 0; j < Nm[i] - 1; ++j)
		{
			k = group[i][j];
			for(unsigned l = 0; l < 3; ++l)
			{
				vrel[k][l] = 2 * (RG.Random() - 0.5 );
				temp[l] += (tempx[j][l] - xcen[i][l]) * vrel[k][l];
			}
		}
     
		for(unsigned l = 0; l < 3; ++l)
		{
			k = Nstart[i]+Nm[i]-1;
			vrel[k][l] = - temp[l] / (tempx[Nm[i]-1][l] - xcen[i][l]);
		}
      
		for(unsigned l = 0; l < 3; ++l)
			sumv[l] = 0;
      
		for(unsigned j = 0; j < Nm[i]; ++j)
		{
			k = Nstart[i] + j;
			for(unsigned l = 0; l < 3; ++l)
				sumv[l] += vrel[k][l] / Nm[i];
		}
      
		for(unsigned j = 0; j < Nm[i]; ++j)
		{
			k = Nstart[i] + j;
			for(unsigned l = 0; l < 3; ++l)
			{
				vrel[k][l] -= sumv[l];
				sumv2[l] += mass[k] * vrel[k][l] * vrel[k][l];
			}
//			std::cout<<k<<"\t"<<vrel[k][0]<<"\t"<<vrel[k][1]<<"\t"<<vrel[k][2]<<"\t\n";
		}
    }

	for(unsigned i = 0; i < 3; ++i)
		sumv2[i] = sqrt((N - 4.0 * Ngroup / 3.0) * Temp / sumv2[i]);

	for(unsigned i = 0; i < N; ++i)
	for(unsigned j = 0; j < 3; ++j)
		vrel[i][j] *= sumv2[j];
  
  
	for(unsigned i = 0; i < Ngroup; ++i)
	for(unsigned j = 0; j < Nm[i]; ++j)
	{
		k = group[i][j];
		for(unsigned l = 0; l < 3; ++l)
			v[k][l] = vcen[i][l] + vrel[k][l];
//		std::cout<<i<<"\t"<<vrel[i][0]<<"\t"<<vrel[i][1]<<"\t"<<vrel[i][2]<<"\t\n";
	}

	
	for(unsigned i = 0; i < 3; ++i)
		sumv[i] = 0;
	
	for(unsigned i = 0; i < Ngroup; ++i)
	for(unsigned j = 0; j < 3; ++j)
		sumv[j] += vcen[i][j];
}



void Atom::center(double* r0, unsigned &igroup)
{
	for(unsigned i = 0; i < 3; ++i)
		r0[i] = 0.0;
	
	double **tempx;
	tempx=new double*[Nm[igroup]];
	
	for(unsigned i = 0; i < Nm[igroup]; ++i)
		tempx[i]=new double[3];
	
	double *ma;
	ma = mass+group[igroup][0];

	for(unsigned i = 0; i < Nm[igroup]; ++i)
	{
		unsigned j = group[igroup][i];
		tempx[i][0] = x[j][0] + boxfold[j][0] * boxx;
		tempx[i][1] = x[j][1] + boxfold[j][1] * boxy;
		tempx[i][2] = x[j][2] + boxfold[j][2] * boxz;
	}

	for(unsigned i = 0; i < Nm[igroup]; ++i)
	for(unsigned j = 0; j < 3; ++j)
		r0[j] += tempx[i][j] * ma[i] / Mass[igroup];
}

double Atom::gyration(unsigned &igroup)
{
	double sumx[3] = {0.0, 0.0, 0.0};
	double **tempx;
	tempx=new double*[Nm[igroup]];
	
	for(unsigned i = 0; i < Nm[igroup]; ++i)
		tempx[i]=new double[3];
	
	for(unsigned i = 0; i < Nm[igroup]; ++i)
	{
		unsigned j = group[igroup][i];

		tempx[i][0] = x[j][0] + boxfold[j][0] * boxx;
		tempx[i][1] = x[j][1] + boxfold[j][1] * boxy;
		tempx[i][2] = x[j][2] + boxfold[j][2] * boxz;	
		
//		std::cout<<boxfold[j][0] <<'\t'<<boxfold[j][1] <<'\t'<<boxfold[j][2] <<"\n";
	}

//	std::cin.get();

	for(unsigned i = 0; i < Nm[igroup]; ++i)
	for(unsigned j = 0; j < 3; ++j)
		sumx[j] += tempx[i][j] / Nm[igroup];
		
	double rg_temp = 0;
	
	for(unsigned i = 0; i < Nm[igroup]; ++i)
	for(unsigned j = 0; j < 3; ++j)
		rg_temp += (tempx[i][j] - sumx[j]) * (tempx[i][j] - sumx[j]);
		


	rg_temp /= Nm[igroup];
	rg_temp = sqrt(rg_temp);


	return rg_temp;
}

void Atom::diffusion(unsigned &loop)
{
	double xtemp, ytemp, ztemp;
	if(loop % dfdt == 0)
	{
		if(loop % dft == 0)
		{
			for(unsigned i = 0; i < Ngroup; ++i)
			for(unsigned j = 0; j < 3; ++j)
				xcen0[i][j] = xcen[i][j];

			stress0 = 0;
			for(unsigned i = 0; i < Ngroup; ++i)
			{
				xtemp = xcen[i][0] - nint(xcen[i][0]/boxx) * boxx;
				ytemp = xcen[i][1] - nint(xcen[i][1]/boxy) * boxy;
				ztemp = xcen[i][2] - nint(xcen[i][2]/boxz) * boxz;

				stress0 += Mass[i] * (xtemp*vcen[i][1] + ytemp*vcen[i][2] + ztemp*vcen[i][0] + 
				ytemp*vcen[i][0] + ztemp*vcen[i][1] + xtemp*vcen[i][2]) / 6.0;
			}
		}
      
		double sumx = 0;
		for(unsigned i = 0; i < Ngroup; ++i)
		for(unsigned j = 0; j < 3; ++j)
			sumx += (xcen[i][j] - xcen0[i][j]) * (xcen[i][j] - xcen0[i][j]);
      
		unsigned it = loop % dft;
		it = it / dfdt;
		unsigned num = loop / dft + 1;
		diff[it][0] += sumx / 6.0 /Ngroup;
		diff[it][1] = diff[it][0] / num;

		double stress = 0;
		for(unsigned i = 0; i < Ngroup; ++i)
		{
			xtemp = xcen[i][0] - nint(xcen[i][0]/boxx) * boxx;
			ytemp = xcen[i][1] - nint(xcen[i][1]/boxy) * boxy;
			ztemp = xcen[i][2] - nint(xcen[i][2]/boxz) * boxz;

			stress += Mass[i] * (xtemp*vcen[i][1] + ytemp*vcen[i][2] + ztemp*vcen[i][0] +
			ytemp*vcen[i][0] + ztemp*vcen[i][1] + xtemp*vcen[i][2]) / 6.0;
		}

		visc[it][0] += (stress - stress0)*(stress - stress0) / 2.0 / boxx / boxy / boxz;
		visc[it][1] = visc[it][0] / num;
	}
  
	if((10*loop) % dfdt == 0)
	{
		if((10*loop) % dft == 0)
		{
			for(unsigned i = 0; i < Ngroup; ++i)
			for(unsigned j = 0; j < 1; ++j)
				vcen0[i][j] = vcen[i][j];
		}
		
		double sum = 0; 
		for(unsigned i = 0; i < Ngroup; ++i)
		for(unsigned j = 0; j < 1; ++j)
			sum += vcen0[i][j] * vcen[i][j];
		
		sum /= Ngroup;

		int it =  loop % (dft/10); 
		it = it * 10 / dfdt;

		int num = loop * 10 / dft + 1;
		v_corr[it][0] += sum;
		v_corr[it][1] = v_corr[it][0] / num;

	}
}


void Atom::radial_distribution()
{
	double *r_dis_temp;
	r_dis_temp = new double[N_radial];
	for(int i = 0; i < N_radial; ++i)
		r_dis_temp[i] = 0;
  
	double vb;
	double rbin = (r_radial / N_radial);
	rbin = rbin * rbin * rbin;
	double rho = Ngroup / boxx / boxy / boxz;

	double r;
	for(unsigned i = 0; i < Ngroup; ++i)
		particle_number[i] = 0;

	for(unsigned i = 0; i < Ngroup -1; ++i)
	for(unsigned j = i + 1; j < Ngroup; ++j)
	{
		double xij = xcen[i][0] - xcen[j][0];
		double yij = xcen[i][1] - xcen[j][1];
		double zij = xcen[i][2] - xcen[j][2];

		xij = xij - nint(xij/boxx) * boxx;
		yij = yij - nint(yij/boxy) * boxy;
		zij = zij - nint(zij/boxz) * boxz;
	
		r = sqrt(xij*xij + yij*yij + zij*zij);
		if(r < boxx / 2.0)
		{
			int index = static_cast<int>(r*N_radial/r_radial);
			if(index < N_radial)
				r_dis_temp[index] += 2.0;
		}
		if(r < 5.5)
		{
			particle_number[i] += 1;
			particle_number[j] += 1;
		}
	}
  
	for(int i = 0; i < N_radial; ++i)
    {
		vb = (4.0/3.0) * PI * ((i+1) * (i+1) * (i+1) - i * i * i) * rbin * rho;
		r_dis_temp[i] = r_dis_temp[i] / (Ngroup * vb);
		r_dis[i] += r_dis_temp[i];
	}
		
	for(unsigned i = 0; i < Ngroup; ++i)
	{
		unsigned index = particle_number[i];
		pnumber_dis[index] += 1;
	}

	delete r_dis_temp;

	++r_dis_num;
}


void Atom::velocity_distribution()
{
	double xtemp;

	for(unsigned i = 0; i < Ngroup; ++i)
	{
		xtemp = xcen[i][0] - nint(xcen[i][0]/boxx) * boxx;
		int index = static_cast<int>((xtemp + 0.5 * boxx) * velocitybin / boxx);
		veldis[index][0] += vcen[i][2];
		veldis[index][1] += 1.0;
	}
}

